function [mu sigma] = computeMLEGaussianUncertainty(data, w, covar_type, reg)


[N D] = size(data.Y);

Y = data.Y;
YYt = data.YYt;
C = data.C;

if sum(w) > 0
    mu = w'*Y/sum(w);
    sigma = (reshape(reshape(YYt+C, D*D, N)*w, D, D) - sum(w)*(mu'*mu))/sum(w);

    switch covar_type
        case 'full'

        case 'diagonal'
            sigma = diag(diag(sigma));
        case 'spherical'
            sigma = eye(D)*mean(diag(sigma));
        otherwise
            error(['Unknown covariance type ', covar_type]);
    end
else
    mu = zeros(1, D);
    sigma = eye(D)*reg;
end

[U S] = svd(sigma);
if any(diag(S) < reg)
    sigma = sigma + eye(D)*reg;
end
